Dane są informacjami o wartościach odżywczych pozycji w menu sieci
kawiarni Starbucks. Napoje podzielone są na kategorie:
(1) Kawa,
(2) Klasyczne napoje espresso,
(3) Charakterystyczne napoje espresso,
(4) Napoje herbaciane Tazo®,
(5) Mrożone napoje wstrząśnięte,
(6) Smoothies,
(7) Kawa mieszana Frappuccino®,
(8) Kawa mieszana Frappuccino® Light,
(9) Frappuccino® Blended Crème.
Zmienne V1, V2, … , V15 określają kolejno zawartość: kalorii, tłuszczu ogółem (g), tłuszczu trans (g), tłuszczu nasyconego (g), sodu (mg), węglowodanów ogółem (g), cholesterolu (mg), błonniku pokarmowego (g), cukrów (g), białka (g), witaminy A (% DV), witaminy C (% DV), wapnia (% DV), żelaza (% DV) i kofeiny (mg).
Źródło danych: https://www.kaggle.com/datasets/henryshan/starbucks.
Celem projektu jest klasyfikacja bez nadzoru i pod nadzorem danych na podstawie wyselekcjonowanych 15 wartości odżywczych.
starbucks <- read.csv('starbucks_data.csv')
head(starbucks)
## Category V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15
## 1 1 3 0.1 0.0 0.0 0 5 0 0 0 0.3 0.0 0 0.00 0 175
## 2 1 4 0.1 0.0 0.0 0 10 0 0 0 0.5 0.0 0 0.00 0 260
## 3 1 5 0.1 0.0 0.0 0 10 0 0 0 1.0 0.0 0 0.00 0 330
## 4 1 5 0.1 0.0 0.0 0 10 0 0 0 1.0 0.0 0 0.02 0 410
## 5 2 70 0.1 0.1 0.0 5 75 10 0 9 6.0 0.1 0 0.20 0 75
## 6 2 100 3.5 2.0 0.1 15 85 10 0 9 6.0 0.1 0 0.20 0 75
any(is.na(starbucks))
## [1] FALSE
Pierwszym krokiem będzie prezentacja histogramów oraz podstawowych statystyk opisowych dla każdej zmiennej.
Podsumowanie:
apply(starbucks[, -1], 2, summary)
## V1 V2 V3 V4 V5 V6 V7
## Min. 3.0000 0.000000 0.000000 0.00000000 0.000000 0.0000 0.00000
## 1st Qu. 122.5000 0.200000 0.100000 0.00000000 0.000000 80.0000 21.00000
## Median 190.0000 2.500000 0.500000 0.00000000 5.000000 127.5000 36.00000
## Mean 201.0413 3.074312 1.399541 0.03944954 6.674312 136.8578 37.16055
## 3rd Qu. 270.0000 5.000000 2.000000 0.10000000 10.000000 180.0000 52.50000
## Max. 510.0000 15.000000 9.000000 0.30000000 40.000000 340.0000 90.00000
## V8 V9 V10 V11 V12 V13 V14
## Min. 0.0000000 0.00000 0.000000 0.0000000 0.00000000 0.0000000 0.00000000
## 1st Qu. 0.0000000 19.00000 4.000000 0.0400000 0.00000000 0.1000000 0.00000000
## Median 0.0000000 33.00000 6.000000 0.0800000 0.00000000 0.2000000 0.03000000
## Mean 0.8440367 33.94037 7.305046 0.1031193 0.03738532 0.2168349 0.08041284
## 3rd Qu. 1.0000000 44.00000 10.000000 0.1500000 0.00000000 0.3000000 0.10000000
## Max. 8.0000000 84.00000 20.000000 0.5000000 1.00000000 0.6000000 0.50000000
## V15
## Min. 0.00000
## 1st Qu. 51.25000
## Median 75.00000
## Mean 89.93119
## 3rd Qu. 143.75000
## Max. 410.00000
Odchylenie standardowe:
sort(apply(starbucks[, -1], 2, sd), decreasing = TRUE)
## V1 V6 V15 V7 V9 V5
## 102.27230816 80.42698815 64.58923623 20.87723395 19.91320160 8.84163698
## V10 V2 V3 V8 V12 V13
## 4.79903238 3.00748387 1.68767551 1.44435537 0.15070662 0.14558188
## V14 V11 V4
## 0.10797386 0.08214517 0.07313689
Wariancja:
sort(apply(starbucks[, -1], 2, var), decreasing = TRUE)
## V1 V6 V15 V7 V9 V5
## 1.045963e+04 6.468500e+03 4.171769e+03 4.358589e+02 3.965356e+02 7.817454e+01
## V10 V2 V3 V8 V12 V13
## 2.303071e+01 9.044959e+00 2.848249e+00 2.086162e+00 2.271249e-02 2.119408e-02
## V14 V11 V4
## 1.165835e-02 6.747829e-03 5.349004e-03
Skośność:
sort(apply(starbucks[, -1], 2, skewness), decreasing = TRUE)
## V12 V8 V11 V4 V5 V3 V14 V2
## 5.3854396 2.7564703 1.8490660 1.7734325 1.6086669 1.5706864 1.5425005 1.0636256
## V15 V13 V10 V6 V9 V1 V7
## 0.8687843 0.6600035 0.6542657 0.4530970 0.4330222 0.3647753 0.3533467
Wniosek: Największymi wariancjami cechują się zmienne V1, V6 oraz V15, które jednocześnie mają największe średnie. Wszystkie zmienne mają dodatnią skośność.
Histogramy:
for(i in 1:15) hist(starbucks[,i+1], main=colnames(starbucks)[i+1], xlab="wartość", ylab="częstość")
Wniosek: Histogramy wyraźnie pokazują dodatnią skośność zmiennych.
Rozkład normalny przypominają jedynie zmienne V1 i V7 - wykonam dla nich
test Shapiro-Wilka, a dla zmiennych V3-4 i V8-12 wykonam test na
jednomodalność.
shapiro.test(starbucks$V1)
##
## Shapiro-Wilk normality test
##
## data: starbucks$V1
## W = 0.98489, p-value = 0.02008
shapiro.test(starbucks$V7)
##
## Shapiro-Wilk normality test
##
## data: starbucks$V7
## W = 0.97972, p-value = 0.003147
Wniosek: Żadna ze zmiennych nie może mieć rozkładu normalnego.
silverman.test(starbucks$V3, k = 1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.2872873
silverman.test(starbucks$V4, k = 1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.008008008
silverman.test(starbucks$V8, k = 1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.002002002
silverman.test(starbucks$V9, k = 1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.5795796
silverman.test(starbucks$V10, k = 1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0.6706707
silverman.test(starbucks$V11, k = 1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0
silverman.test(starbucks$V12, k = 1)
## Silvermantest: Testing the null hypothesis that the number of modes is <= 1
## The resulting p-value is 0
Wniosek: Przyjmując kryterium 5% błędu otrzymujemy informację, że zmienne V4, V8, V11 i V12 są wielomodalne. Zobaczmy jak wyglądają wykresy wartości p dla tych zmiennych:
silverman.plot(starbucks$V4)
silverman.plot(starbucks$V8)
silverman.plot(starbucks$V11)
silverman.plot(starbucks$V12)
Wniosek: Z wykresów wynika (nadal przyjmując kryterium 5% błędu), że
zmienna V4 powinna być czteromodalna, V8 dwumodalna, V11 dwumodalna a
V12 trójmodalna lub czteromodalna.
Kolejny test - Grubbsa - określi czy próba posiada wartości odstające.
apply(starbucks[, -1], 2, function(x) grubbs.test(x))
## $V1
##
## Grubbs test for one outlier
##
## data: x
## G = 3.02094, U = 0.95775, p-value = 0.2494
## alternative hypothesis: highest value 510 is an outlier
##
##
## $V2
##
## Grubbs test for one outlier
##
## data: x
## G = 3.96534, U = 0.92721, p-value = 0.005931
## alternative hypothesis: highest value 15 is an outlier
##
##
## $V3
##
## Grubbs test for one outlier
##
## data: x
## G = 4.50351, U = 0.90611, p-value = 0.0004397
## alternative hypothesis: highest value 9 is an outlier
##
##
## $V4
##
## Grubbs test for one outlier
##
## data: x
## G = 3.56250, U = 0.94124, p-value = 0.03307
## alternative hypothesis: highest value 0.3 is an outlier
##
##
## $V5
##
## Grubbs test for one outlier
##
## data: x
## G = 3.76918, U = 0.93423, p-value = 0.01402
## alternative hypothesis: highest value 40 is an outlier
##
##
## $V6
##
## Grubbs test for one outlier
##
## data: x
## G = 2.52580, U = 0.97047, p-value = 1
## alternative hypothesis: highest value 340 is an outlier
##
##
## $V7
##
## Grubbs test for one outlier
##
## data: x
## G = 2.53096, U = 0.97034, p-value = 1
## alternative hypothesis: highest value 90 is an outlier
##
##
## $V8
##
## Grubbs test for one outlier
##
## data: x
## G = 4.95443, U = 0.88636, p-value = 3.728e-05
## alternative hypothesis: highest value 8 is an outlier
##
##
## $V9
##
## Grubbs test for one outlier
##
## data: x
## G = 2.51389, U = 0.97074, p-value = 1
## alternative hypothesis: highest value 84 is an outlier
##
##
## $V10
##
## Grubbs test for one outlier
##
## data: x
## G = 2.6453, U = 0.9676, p-value = 0.8416
## alternative hypothesis: highest value 20 is an outlier
##
##
## $V11
##
## Grubbs test for one outlier
##
## data: x
## G = 4.83146, U = 0.89193, p-value = 7.508e-05
## alternative hypothesis: highest value 0.5 is an outlier
##
##
## $V12
##
## Grubbs test for one outlier
##
## data: x
## G = 6.38734, U = 0.81112, p-value = 2.027e-09
## alternative hypothesis: highest value 1 is an outlier
##
##
## $V13
##
## Grubbs test for one outlier
##
## data: x
## G = 2.63196, U = 0.96793, p-value = 0.8765
## alternative hypothesis: highest value 0.6 is an outlier
##
##
## $V14
##
## Grubbs test for one outlier
##
## data: x
## G = 3.88601, U = 0.93009, p-value = 0.008446
## alternative hypothesis: highest value 0.5 is an outlier
##
##
## $V15
##
## Grubbs test for one outlier
##
## data: x
## G = 4.95545, U = 0.88631, p-value = 3.706e-05
## alternative hypothesis: highest value 410 is an outlier
Wniosek: Przyjmując kryterium 5% błędu mamy, że jedynie zmienne V1, V6, V7, V9 i V10 nie mają wartości odstających.
Dokonujemy obliczeń macierzy korelacji Pearsona oraz współczynnika zależności Kendalla, aby uniezależnić się od założenia normalności rozkładu danych.
heatmaply_cor(cor(starbucks[, -1]))
Wniosek: Mapa cieplna jest bardzo zróżnicowana. Najsilniej skorelowane są zbiory zmiennych {V1, V6, V7, V9}, {V4, V5}, {V10, V11, V13} oraz {V8, V12}. Słabiej skorelowane są zbiory takie jak {V3, V4, V5} i {V2, V3, V5, V10, V13}. Najsłabiej są skorelowane np. zmienne V7 i V14.
heatmaply_cor(cor(starbucks[, -1], method = 'kendall'))
Wniosek: Macierz korelacji Kendalla jest podobna do macierzy Pearsona, widać na niej jednak nieco słabsze skorelowanie między zmiennymi.\
prcomp(starbucks[, -1]) -> pca.starbucks
summary(pca.starbucks)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 126.1494 65.4197 40.69682 11.77600 6.70131 3.09895
## Proportion of Variance 0.7218 0.1941 0.07512 0.00629 0.00204 0.00044
## Cumulative Proportion 0.7218 0.9159 0.99102 0.99731 0.99935 0.99979
## PC7 PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 1.89625 0.84376 0.46864 0.42542 0.07908 0.05196 0.03425
## Proportion of Variance 0.00016 0.00003 0.00001 0.00001 0.00000 0.00000 0.00000
## Cumulative Proportion 0.99995 0.99998 0.99999 1.00000 1.00000 1.00000 1.00000
## PC14 PC15
## Standard deviation 0.02352 0.01781
## Proportion of Variance 0.00000 0.00000
## Cumulative Proportion 1.00000 1.00000
prcomp(starbucks[, -1], scale. = TRUE) -> pca.starbucks.sc
summary(pca.starbucks.sc)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.4715 1.5835 1.4521 1.16984 1.07465 0.88856 0.54014
## Proportion of Variance 0.4072 0.1672 0.1406 0.09124 0.07699 0.05264 0.01945
## Cumulative Proportion 0.4072 0.5744 0.7150 0.80620 0.88319 0.93583 0.95528
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.51786 0.43590 0.33292 0.24555 0.15849 0.11681 0.04168
## Proportion of Variance 0.01788 0.01267 0.00739 0.00402 0.00167 0.00091 0.00012
## Cumulative Proportion 0.97315 0.98582 0.99321 0.99723 0.99890 0.99981 0.99993
## PC15
## Standard deviation 0.03235
## Proportion of Variance 0.00007
## Cumulative Proportion 1.00000
df.pca <- data.frame(pc1 = pca.starbucks$x[, 1], pc2 = pca.starbucks$x[, 2], pc3 = pca.starbucks$x[, 3], kl = as.factor(starbucks$Category))
plot(x = df.pca[, 1], y = df.pca[, 2], col = as.factor(starbucks[, 1]))
plot_ly(df.pca, x=~pc1, y=~pc2, z=~pc3, color=~kl, type='scatter3d', colors = viridis(9))
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
pca.starbucks$rotation[, 1] -> pca.rot1
sort(pca.rot1)
## V15 V4 V12 V11 V14
## -7.695161e-05 1.752291e-04 1.761441e-04 2.217498e-04 3.216063e-04
## V13 V8 V3 V2 V10
## 5.175786e-04 3.402807e-03 8.268791e-03 1.350765e-02 1.877526e-02
## V5 V9 V7 V6 V1
## 2.470073e-02 1.443410e-01 1.536808e-01 5.770364e-01 7.882585e-01
pca.starbucks$rotation[, 2] -> pca.rot2
sort(pca.rot2)
## V15 V6 V5 V2 V3
## -9.796140e-01 -1.628696e-01 -5.026304e-03 -3.768940e-03 -2.980319e-03
## V13 V4 V14 V11 V10
## -2.306455e-04 -5.301965e-06 9.095648e-06 2.818636e-05 3.950922e-04
## V12 V8 V9 V7 V1
## 6.530707e-04 4.909092e-03 3.095668e-02 3.964840e-02 1.059556e-01
pca.starbucks$rotation[, 3] -> pca.rot3
sort(pca.rot3)
## V6 V4 V11 V14 V13
## -0.7984239836 0.0002377182 0.0002803016 0.0004738972 0.0007056207
## V12 V3 V8 V2 V10
## 0.0007911645 0.0073478188 0.0094140046 0.0211193312 0.0300239710
## V5 V9 V7 V15 V1
## 0.0322583642 0.0425003221 0.0735806015 0.1974442763 0.5601627919
Wniosek: Dla niewyskalowanego PCA mamy, że już pierwsze dwie składowe zawierają 92% wariancji, a trzy składowe 99%, natomiast dla wyskalowanego PCA trzy składowe zawierają zaledwie 72%. Dalej dla niewyskalowanego PCA: PC1 i PC2 zdominowane są przez zmienną V15. Niestety na wykresach nie widać aby ta zmienna oddzielała którykolwiek klaster od reszty.
set.seed(10)
Rtsne(starbucks[, -1], theta = 0, dims = 2) -> starbucks.tsne2
summary(starbucks.tsne2$Y)
## V1 V2
## Min. :-8.3469 Min. :-16.0396
## 1st Qu.:-5.0850 1st Qu.: -6.5213
## Median : 0.3785 Median : 0.1484
## Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 4.1747 3rd Qu.: 6.9102
## Max. : 9.8827 Max. : 14.1478
sd(starbucks.tsne2$Y)
## [1] 6.805308
tsne2.y.df <- as.data.frame(starbucks.tsne2$Y)
ggplot(tsne2.y.df, aes(x = V1, y = V2, col = starbucks$Category)) +
geom_point() + labs(x = "t-SNE 1", y = "t-SNE 2", col = "Category") + theme_classic()
Rtsne(starbucks[, -1], theta = 0, dims = 3) -> starbucks.tsne3
summary(starbucks.tsne3$Y)
## V1 V2 V3
## Min. :-17.132 Min. :-25.8961 Min. :-22.5534
## 1st Qu.: -5.677 1st Qu.:-10.0151 1st Qu.: -8.3242
## Median : -0.433 Median : -0.9437 Median : 0.5223
## Mean : 0.000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 8.073 3rd Qu.: 10.4523 3rd Qu.: 8.4342
## Max. : 13.631 Max. : 26.3293 Max. : 17.9610
sd(starbucks.tsne3$Y)
## [1] 11.15663
tsne3.y.df <- as.data.frame(starbucks.tsne3$Y)
plot_ly(tsne3.y.df, x = ~V1, y = ~V2, z = ~V3, color = starbucks$Category, type = 'scatter3d', colors = "Set3")
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Wniosek: T-SNE jedynie trochę lepiej radzi sobie z odróżnianiem kategorii, nadal jest to raczej chaotyczna mieszanka.
Sprawdzę jak działa k-średnich dla ilości klastrów równej ilości kategorii kawy.
starbucks.km <- kmeans(starbucks[, -1], centers = 9)
df.pca2 <- data.frame(PC1 = pca.starbucks$x[, 1], PC2 = pca.starbucks$x[, 2], PC3 = pca.starbucks$x[, 3], kl = as.factor(starbucks.km$cluster))
plot_ly(df.pca2, x = ~PC1, y = ~PC2, z = ~PC3, color = starbucks.km$cluster, type = 'scatter3d')
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Wniosek: Wizualnie k-średnich ładnie podzieliło wykres.
Porównam jednak w tabeli czy jest to jakkolwiek wierne faktycznemu podziałowi.
table(starbucks$Category, starbucks.km$cluster)
##
## 1 2 3 4 5 6 7 8 9
## 1 0 0 3 1 0 0 0 0 0
## 2 17 12 2 8 7 12 0 0 0
## 3 7 1 0 0 7 4 9 6 6
## 4 13 0 0 0 8 3 6 3 3
## 5 1 0 1 7 3 0 0 0 0
## 6 0 0 0 0 0 0 8 0 0
## 7 0 1 0 0 7 0 1 14 13
## 8 1 7 0 0 0 0 0 3 1
## 9 0 0 0 0 0 0 9 3 0
Wniosek: Jak widać z tabeli k-średnich słabo poradziło sobie z klasteryzacją, np. zmienna V2 została w głównej mierze podzielona między kl1 a kl5, podobnie zmienna V7 pomiędzy kl2 i kl7.
Użyję metody łokcia do znalezienia optymalnej ilości klastrów.
wss <- rep(NA, 9)
for(i in 1:9) wss[i] <- kmeans(starbucks[, -1], centers = i + 1)$tot.withinss
plot(1:9, wss)
Wniosek: Według wykresu optymalna ilość klastrów to 5.
Zatem teraz sprawdzę k-średnich dla ilości klastrów równej 5.
starbucks.km5 <- kmeans(starbucks[, -1], centers = 5)
df.pca3 <- data.frame(PC1 = pca.starbucks$x[, 1], PC2 = pca.starbucks$x[, 2], PC3 = pca.starbucks$x[, 3], kl = as.factor(starbucks.km$cluster))
plot_ly(df.pca3, x = ~PC1, y = ~PC2, z = ~PC3, color = starbucks.km5$cluster,type = 'scatter3d')
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
table(starbucks$Category, starbucks.km5$cluster)
##
## 1 2 3 4 5
## 1 0 0 0 0 4
## 2 29 0 2 21 6
## 3 11 11 8 10 0
## 4 17 8 4 7 0
## 5 8 0 0 0 4
## 6 0 8 0 0 0
## 7 0 1 16 19 0
## 8 3 0 1 8 0
## 9 0 10 2 0 0
Wniosek: Nadal jest to bardzo nieprecyzyjny podział.
plot(agnes(dist(starbucks[, -1], method = 'minkowski', p = 1), diss = T))
plot(diana(starbucks[, -1]))
Wniosek: Metoda aglomeracyjna i podziałowa dają względem siebie bardzo
zbliżone rezultaty.
rfcv(trainx = starbucks[, -1], trainy = as.factor(starbucks$Category)) -> rf.starbucks
rf.starbucks
## $n.var
## [1] 15 8 4 1
##
## $error.cv
## 15 8 4 1
## 0.1513761 0.1238532 0.1376147 0.4357798
##
## $predicted
## $predicted$`15`
## [1] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 2 2 2 2 3 2 2 3 2 2 3 3 2 3 2
## [38] 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2
## [75] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4
## [112] 4 3 4 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5
## [149] 5 5 6 6 6 6 6 6 6 6 9 7 8 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [186] 7 7 7 7 7 7 7 7 7 8 8 8 8 8 7 8 8 8 8 7 7 9 9 9 9 9 9 9 9 7 7 7 7
## Levels: 1 2 3 4 5 6 7 8 9
##
## $predicted$`8`
## [1] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 2 3 2 3 2 2 2
## [38] 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 3 2 2 2 3
## [75] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 3 4 4 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4
## [112] 4 3 4 4 4 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 4 4 4 3 4 3 4 5 5 5 5 5 5 5 5 5 5
## [149] 5 5 6 6 6 6 6 6 6 6 9 7 8 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [186] 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 7 7 9 9 9 9 9 9 9 9 9 9 7 9
## Levels: 1 2 3 4 5 6 7 8 9
##
## $predicted$`4`
## [1] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 2 3 2 2 2
## [38] 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 3 2 3 3 2 3 2 3
## [75] 3 3 2 7 3 7 3 3 2 3 3 3 4 3 4 3 3 4 3 3 3 3 3 3 4 3 3 3 4 4 4 4 4 4 4 4 4
## [112] 4 4 4 4 4 4 4 4 4 4 4 4 4 7 4 4 4 4 4 4 4 4 4 4 4 9 4 5 5 5 5 5 5 5 5 5 5
## [149] 5 5 6 6 6 6 6 6 6 6 9 7 7 7 7 7 7 7 7 4 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [186] 3 3 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 3 9 9 7 9 9 9 9 9 9 9 7 9
## Levels: 1 2 3 4 5 6 7 8 9
##
## $predicted$`1`
## [1] 2 5 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 4 7 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 5 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 4 3 4 3 3 3 4 9 4 4 4 4 4 7 7 7 2 2 7
## [112] 8 8 8 3 4 3 4 4 4 4 4 4 7 7 7 4 9 9 9 9 4 4 9 9 9 9 9 8 8 2 5 5 5 5 5 5 5
## [149] 5 5 6 6 6 4 4 9 9 4 4 7 7 7 4 2 4 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [186] 2 2 2 7 7 7 7 7 7 7 7 4 7 4 7 7 5 4 7 7 5 9 9 9 4 4 4 4 9 4 4 9 4
## Levels: 1 2 3 4 5 6 7 8 9
Wniosek: Błąd walidacji krzyżowej jest najniższy przy 8 cechach - wybierzemy te 8 najlepszych cech.
randomForest(starbucks[, -1], as.factor(starbucks$Category), importance = T) -> rf.starbucks.imp
varImpPlot(rf.starbucks.imp)
Wniosek: Najlepsze są cechy o numerach 15, 9, 7, 6, 10, 13, 1 i 11.
summary(as.factor(starbucks$Category))
## 1 2 3 4 5 6 7 8 9
## 4 58 40 36 12 8 36 12 12
model <- rpart(as.factor(starbucks$Category) ~ ., method = "class", data=starbucks[, -1])
printcp(model)
##
## Classification tree:
## rpart(formula = as.factor(starbucks$Category) ~ ., data = starbucks[,
## -1], method = "class")
##
## Variables actually used in tree construction:
## [1] V12 V13 V15 V3 V6 V7 V8 V9
##
## Root node error: 160/218 = 0.73394
##
## n= 218
##
## CP nsplit rel error xerror xstd
## 1 0.11875 0 1.00000 1.00000 0.040778
## 2 0.08750 2 0.76250 0.77500 0.045701
## 3 0.07500 3 0.67500 0.73750 0.045983
## 4 0.06875 4 0.60000 0.73125 0.046016
## 5 0.05625 5 0.53125 0.64375 0.046070
## 6 0.05000 6 0.47500 0.61250 0.045905
## 7 0.03750 8 0.37500 0.60000 0.045811
## 8 0.02500 9 0.33750 0.50000 0.044477
## 9 0.01000 11 0.28750 0.47500 0.043975
fancyRpartPlot(model, )
table(predict(model, starbucks, type="class"), starbucks$Category)
##
## 1 2 3 4 5 6 7 8 9
## 1 0 0 0 0 0 0 0 0 0
## 2 4 47 2 9 0 0 0 0 0
## 3 0 7 29 4 0 0 0 0 0
## 4 0 4 8 23 0 0 0 0 0
## 5 0 0 0 0 12 0 0 0 0
## 6 0 0 1 0 0 6 0 0 0
## 7 0 0 0 0 0 2 35 4 0
## 8 0 0 0 0 0 0 0 8 0
## 9 0 0 0 0 0 0 1 0 12
set.seed(10)
randomForest(starbucks[, c(16, 10, 8, 7, 11, 14, 2, 12)], as.factor(starbucks$Category)) -> rf.starbucks.sel
rf.starbucks.sel$confusion
## 1 2 3 4 5 6 7 8 9 class.error
## 1 3 1 0 0 0 0 0 0 0 0.25000000
## 2 0 55 3 0 0 0 0 0 0 0.05172414
## 3 0 8 27 4 1 0 0 0 0 0.32500000
## 4 0 0 5 31 0 0 0 0 0 0.13888889
## 5 0 0 0 0 12 0 0 0 0 0.00000000
## 6 0 0 0 0 0 8 0 0 0 0.00000000
## 7 0 0 0 0 0 0 34 1 1 0.05555556
## 8 0 0 0 0 0 0 2 10 0 0.16666667
## 9 0 0 0 0 0 0 0 0 12 0.00000000
1 - mean(rf.starbucks.sel$confusion[, 10])
## [1] 0.8902405
Wniosek: Metodą lasów losowych osiągnęliśmy średnią dokładność na poziomie 89%, co stanowi ulepszenie w porównaniu z drzewami decyzyjnymi. Największy błąd występuje w klasie V3 i wynosi 33%.
party <- createDataPartition(starbucks[, 1], p = 0.7, list = F)
starbucks.test <- starbucks[-party,]
starbucks.train <- starbucks[party,]
prop.table(table(starbucks$Category))
##
## 1 2 3 4 5 6 7
## 0.01834862 0.26605505 0.18348624 0.16513761 0.05504587 0.03669725 0.16513761
## 8 9
## 0.05504587 0.05504587
prop.table(table(starbucks.train$Category))
##
## 1 2 3 4 5 6 7
## 0.01935484 0.26451613 0.19354839 0.15483871 0.05806452 0.03225806 0.16774194
## 8 9
## 0.05806452 0.05161290
preparty <- preProcess(starbucks.train[, -1], method=c("center", "scale"))
starbucks.train.sc <- predict(preparty, starbucks.train)
starbucks.test.sc <- predict(preparty, starbucks.test)
apply(starbucks.train.sc[,-1],2,sd)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
summary(starbucks.train.sc)
## Category V1 V2 V3
## Min. :1.000 Min. :-1.95136 Min. :-1.0067 Min. :-0.8251
## 1st Qu.:2.000 1st Qu.:-0.74160 1st Qu.:-0.9407 1st Qu.:-0.7665
## Median :4.000 Median :-0.09705 Median :-0.1816 Median :-0.5322
## Mean :4.335 Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.:7.000 3rd Qu.: 0.69624 3rd Qu.: 0.5609 3rd Qu.: 0.3466
## Max. :9.000 Max. : 3.07610 Max. : 3.9438 Max. : 4.4474
## V4 V5 V6 V7
## Min. :-0.5387 Min. :-0.7773 Min. :-1.6748 Min. :-1.76838
## 1st Qu.:-0.5387 1st Qu.:-0.7773 1st Qu.:-0.7316 1st Qu.:-0.71544
## Median :-0.5387 Median :-0.1867 Median :-0.1535 Median :-0.04539
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.9010 3rd Qu.: 0.4039 3rd Qu.: 0.5767 3rd Qu.: 0.64859
## Max. : 3.7804 Max. : 3.3569 Max. : 2.4631 Max. : 2.53909
## V8 V9 V10 V11
## Min. :-0.5670 Min. :-1.69731 Min. :-1.5470 Min. :-1.2910
## 1st Qu.:-0.5670 1st Qu.:-0.69279 1st Qu.:-0.6898 1st Qu.:-0.7847
## Median :-0.5670 Median :-0.09008 Median :-0.2613 Median :-0.2785
## Mean : 0.0000 Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.1654 3rd Qu.: 0.51263 3rd Qu.: 0.5958 3rd Qu.: 0.6075
## Max. : 4.5598 Max. : 2.52166 Max. : 2.5244 Max. : 5.0375
## V12 V13 V14 V15
## Min. :-0.2256 Min. :-1.5027 Min. :-0.7346 Min. :-1.4288
## 1st Qu.:-0.2256 1st Qu.:-0.8045 1st Qu.:-0.7346 1st Qu.:-0.6249
## Median :-0.2256 Median :-0.1063 Median :-0.5555 Median :-0.2230
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.:-0.2256 3rd Qu.: 0.5919 3rd Qu.: 0.1613 3rd Qu.: 0.9024
## Max. : 6.6581 Max. : 2.6864 Max. : 3.7449 Max. : 3.8766
Wniosek: Podzieliłam dane na część treningową (70%) i testową (30%), tak że ich procentowy podział klas pozostał bez większych zmian. Następnie wyskalowałam/wystandaryzowałam je tak, że w danych treningowych odchylenie standardowe wynosi 1, a średnia 0.
Wykorzystam ten podział danych do KNN (k najbliższych sąsiadów), policzę również optymalne k dzięki pętli for.
set.seed(10)
długość <- length(starbucks.train.sc$Category)
max.dokładność <- 0
max.k <- 0
for (i in c(2:długość)){
prognoza <- knn(starbucks.train.sc[, -1],starbucks.test.sc[, -1],starbucks.train.sc$Category,k = i)
dokładność <- mean(prognoza == starbucks.test.sc$Category)
if (max.dokładność < dokładność){
max.k = i
max.dokładność = dokładność
}
}
c(max.k, max.dokładność)
## [1] 3.0000000 0.8095238
knn(starbucks.train.sc[, -1],starbucks.test.sc[, -1], starbucks.train.sc$Category, k = max.k) -> starbucks.knn
conf_matrix.knn <- table(starbucks.knn, starbucks.test.sc$Category)
conf_matrix.knn
##
## starbucks.knn 1 2 3 4 5 6 7 8 9
## 1 1 1 0 0 0 0 0 0 0
## 2 0 12 4 1 0 0 0 0 0
## 3 0 4 5 4 0 0 0 0 0
## 4 0 0 1 7 0 0 0 0 0
## 5 0 0 0 0 3 0 0 0 0
## 6 0 0 0 0 0 3 0 0 0
## 7 0 0 0 0 0 0 10 0 1
## 8 0 0 0 0 0 0 0 3 0
## 9 0 0 0 0 0 0 0 0 3
mean(starbucks.knn == starbucks.test.sc$Category)
## [1] 0.7460317
Wniosek: Najlepsze k dla moich danych to 2, a na macierzy widać że knn średnio sobie poradziło.
Wykonam teraz Naiwny klasyfikator Bayes’a, również przy użyciu stworzonych wcześniej danych treningowych i testowych.
model.nb <- naiveBayes(Category~., data = starbucks.train.sc)
model.nb
##
## Naive Bayes Classifier for Discrete Predictors
##
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
##
## A-priori probabilities:
## Y
## 1 2 3 4 5 6 7
## 0.01935484 0.26451613 0.19354839 0.15483871 0.05806452 0.03225806 0.16774194
## 8 9
## 0.05806452 0.05161290
##
## Conditional probabilities:
## V1
## Y [,1] [,2]
## 1 -1.9414457 0.009916107
## 2 -0.5771346 0.793167598
## 3 0.5078326 1.076330422
## 4 -0.1301036 0.739552706
## 5 -0.8242311 0.451700157
## 6 0.7953997 0.140234926
## 7 0.7763303 0.789819849
## 8 -0.2182468 0.577730830
## 9 0.3243846 0.692608034
##
## V2
## Y [,1] [,2]
## 1 -0.97370076 0.0000000
## 2 0.05181320 1.0076355
## 3 0.72046333 1.2669437
## 4 -0.12524358 0.7538134
## 5 -0.77568158 0.2640256
## 6 -0.34664003 0.4811001
## 7 -0.05341932 0.8134777
## 8 -0.51165601 0.6475695
## 9 -0.38789402 0.7400268
##
## V3
## Y [,1] [,2]
## 1 -0.8250782 0.0000000
## 2 -0.1449416 0.8320957
## 3 0.7664320 1.3069830
## 4 -0.2905065 0.6331589
## 5 -0.7014026 0.1974620
## 6 -0.4032792 0.3055326
## 7 0.1843554 1.0532350
## 8 -0.1806630 0.8884582
## 9 -0.3051523 0.6653312
##
## V4
## Y [,1] [,2]
## 1 -0.53873064 0.0000000
## 2 0.05822278 1.1140616
## 3 0.18112496 1.2396401
## 4 0.12113699 1.0381182
## 5 -0.37876273 0.4799037
## 6 -0.53873064 0.0000000
## 7 0.01500443 0.8222173
## 8 -0.53873064 0.0000000
## 9 0.18112496 1.0883194
##
## V5
## Y [,1] [,2]
## 1 -0.7773001 0.0000000
## 2 0.1157958 1.1215251
## 3 0.4038912 1.2789710
## 4 0.1332015 1.0732825
## 5 -0.5804349 0.4176142
## 6 -0.1867044 0.0000000
## 7 -0.2321349 0.6233082
## 8 -0.5148131 0.3112712
## 9 -0.1128800 0.6650058
##
## V6
## Y [,1] [,2]
## 1 -1.5733828 0.03513289
## 2 -0.3891940 0.62695087
## 3 0.0290519 0.97821058
## 4 -0.5388996 0.46157812
## 5 -1.3042820 0.23068517
## 6 -0.1048224 0.14528204
## 7 1.1337489 0.73016100
## 8 1.1649550 0.70441227
## 9 0.6299649 0.65574860
##
## V7
## Y [,1] [,2]
## 1 -1.76838012 0.0000000
## 2 -0.79949034 0.6654875
## 3 0.29442124 0.9239913
## 4 -0.15906003 0.6112599
## 5 -0.57185961 0.4068170
## 6 0.85439284 0.1490617
## 7 1.10142830 0.7963305
## 8 0.02374171 0.5340885
## 9 0.63662611 0.6767333
##
## V8
## Y [,1] [,2]
## 1 -0.56701657 0.0000000
## 2 -0.06684342 0.7220519
## 3 -0.02992587 0.6916995
## 4 -0.23133489 0.5705570
## 5 -0.56701657 0.0000000
## 6 4.26679970 0.4011500
## 7 -0.03180381 0.6059957
## 8 -0.07875230 0.6342739
## 9 -0.47546702 0.2589412
##
## V9
## Y [,1] [,2]
## 1 -1.69730775 0.0000000
## 2 -0.85694449 0.5996699
## 3 0.25982396 0.9215496
## 4 -0.08589691 0.6177917
## 5 -0.45282422 0.4222984
## 6 0.18113689 0.2058645
## 7 1.18874356 0.7714827
## 8 0.07175628 0.4763366
## 9 0.76375608 0.7072493
##
## V10
## Y [,1] [,2]
## 1 -1.41840270 0.07726015
## 2 0.32929493 0.91010271
## 3 0.59583973 1.13944092
## 4 -0.02914684 0.65030504
## 5 -1.08983833 0.30379590
## 6 2.01009505 0.28748828
## 7 -0.63215586 0.31816346
## 8 -0.57080186 0.30513643
## 9 -0.63627665 0.27465264
##
## V11
## Y [,1] [,2]
## 1 -1.29102775 0.0000000
## 2 0.32661582 0.8911882
## 3 0.39659023 0.8802654
## 4 0.08543566 0.6489549
## 5 -0.89725023 0.2861496
## 6 1.59479899 3.1444766
## 7 -0.59001721 0.2801336
## 8 -0.47534573 0.2109522
## 9 -0.46831399 0.2620276
##
## V12
## Y [,1] [,2]
## 1 -0.22560925 0.00000000
## 2 -0.21217755 0.04135841
## 3 -0.17053928 0.07754429
## 4 -0.05064735 0.36310884
## 5 -0.22560925 0.00000000
## 6 4.93720085 2.40931138
## 7 -0.22560925 0.00000000
## 8 -0.22560925 0.00000000
## 9 0.11857809 0.22077105
##
## V13
## Y [,1] [,2]
## 1 -1.5026628 0.0000000
## 2 0.4624575 0.9513657
## 3 0.7082394 1.1680382
## 4 0.1846047 0.7611138
## 5 -1.0372097 0.2792719
## 6 -0.6648472 0.3122354
## 7 -0.6594766 0.2792182
## 8 -0.6648472 0.2035526
## 9 -0.5426658 0.3094351
##
## V14
## Y [,1] [,2]
## 1 -0.73463724 0.0000000
## 2 0.05200585 1.1435463
## 3 0.19709780 1.1580197
## 4 -0.30535227 0.5751615
## 5 -0.61518403 0.2003291
## 6 0.55545742 0.9258587
## 7 0.24396021 1.0245962
## 8 0.20107957 1.1750555
## 9 -0.48826499 0.1900489
##
## V15
## Y [,1] [,2]
## 1 2.6708264 1.2478967
## 2 0.5553523 0.7936152
## 3 -0.2390779 0.9132074
## 4 -0.7655937 0.6418702
## 5 0.6880187 0.5131377
## 6 -1.3323015 0.1320845
## 7 0.1572776 0.5354125
## 8 0.2860982 0.5019985
## 9 -1.4287624 0.0000000
pred.nb <- predict(model.nb, starbucks.test.sc)
pred.nb
## [1] 1 2 2 8 8 2 3 5 2 3 3 5 5 2 2 5 2 2 2 2 2 4 3 3 8 3 3 9 8 3 3 4 4 4 4 4 4 4
## [39] 5 8 1 1 5 6 6 6 8 7 7 8 8 8 8 7 7 7 8 8 8 9 9 9 8
## Levels: 1 2 3 4 5 6 7 8 9
conf_matrix.bayes <- table(pred.nb, starbucks.test.sc$Category)
conf_matrix.bayes
##
## pred.nb 1 2 3 4 5 6 7 8 9
## 1 1 0 0 0 2 0 0 0 0
## 2 0 8 3 0 0 0 0 0 0
## 3 0 3 4 2 0 0 0 0 0
## 4 0 0 1 7 0 0 0 0 0
## 5 0 4 0 1 1 0 0 0 0
## 6 0 0 0 0 0 3 0 0 0
## 7 0 0 0 0 0 0 5 0 0
## 8 0 2 1 2 0 0 5 3 1
## 9 0 0 1 0 0 0 0 0 3
mean(pred.nb == starbucks.test.sc$Category)
## [1] 0.5555556
Wniosek: Naiwny klasyfikator Bayes’a również poradził sobie jedynie w około 60%.
Celem projektu było przeprowadzenie klasyfikacji danych zarówno bez nadzoru, jak i pod nadzorem, opierając się na 15 wartościach odżywczych.
Przeprowadziłam analizę rozkładów poszczególnych zmiennych,
identyfikując wariancję, skośność i wielomodalność:
- wszystkie zmienne wykazują dodatnią skośność oraz nie są podobne do
rozkładu normalnego,
- zmienne V4, V8, V11 i V12 są prawdopodobnie wielomodalne,
- 2/3 zmiennych posiada wartości odstające.
Zbadałam macierz korelacji Pearsona i współczynnika zależności
Kendalla oraz zastosowałam PCA i t-SNE do redukcji wymiaru danych:
- macierze wykazały, że spora część zmiennych wykazuje silną
korelację,
- metody redukcji wymiaru niestety nie zapewniły jednoznacznej separacji
danych.
Przeprowadziłam analizę skupień za pomocą k-średnich, aglomeracyjnej i podziałowej, z identyfikacją optymalnej liczby klastrów. Wyniki analizy skupień były jednak niezadowalające, co sugeruje trudności w jednoznacznym podziale danych na grupy.
Zastosowałam algorytm lasów losowych do klasyfikacji, wykorzystując istotność cech do wyboru najlepszych predyktorów. Uzyskałam dokładność klasyfikacji na poziomie 89%.
Przeprowadziłam klasyfikację za pomocą k najbliższych sąsiadów (KNN) i naiwnego klasyfikatora Bayes’a. Oba modele osiągnęły dokładność na poziomie około 60%.
Pomimo zastosowania różnorodnych metod analizy danych i klasyfikacji, uzyskane rezultaty sugerują trudności w jednoznacznym modelowaniu i klasyfikacji danych dotyczących wartości odżywczych w menu Starbucks.